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Abstract 

In this paper we perform an analytical and numerical study of Extreme Value distributions 
in discrete dynamical systems that have a singular measure. Using the block maxima approach 



described in Faranda et al. 2011 we show that, numerically, the Extreme Value distribution for 
these maps can be associated to the Generalised Extreme Value family where the parameters 
scale with the information dimension. The numerical analysis are performed on a few low 
dimensional maps. For the middle third Cantor set and the Sierpinskij triangle obtained using 
Iterated Function Systems, experimental parameters show a very good agreement with the 
theoretical values. For strange attractors like Lozi and Henon maps a slower convergence to 
the Generalised Extreme Value distribution is observed. Even in presence of large statistics the 
observed convergence is slower if compared with the maps which have an absolute continuous 
invariant measure. Nevertheless and within the uncertainty computed range, the results are in 
good agreement with the theoretical estimates. 

The existence of extreme value laws for dynamical systems preserving an absolutely 
continuous invariant measure or a singular continuous invariant measure has been re- 
cently proven if strong mixing properties or exponential hitting time statistics on balls 
are satisfied. In our previous work we have shown that there exists an algorithmic way 
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to study extrema by using a block-maxima approach for dynamical systems which pos- 
sess an absolutely continuous invariant measure and satisfy certain mixing properties. 
In this work we test our algorithm for maps that do not have an absolutely contin- 
uous invariant measure showing that the cumulative distribution function of maxima 
is related to the scaling of the measure of a ball centered around generic points. The 
scaling exponent turns out to be the Hausdorff dimension of the measure (also known 
as information dimension). Even if we cannot estimate analytically the asymptotic 
behavior of the measure of the balls, the agreement with the numerical simulations 
we have carried out for different maps suggests the validity of our proposed scaling in 
terms of the information dimension. Our conjecture has been tested with numerical 
experiments on different low dimensional maps such as the middle third Cantor set, 
the Sierpinskij triangle, Iterated Function System (IFS) with non-uniform weights and 
strange attractors such as Lozi and Henon. In all cases considered, there is a good 
agreement between the theoretical parameters and the experimental ones although, 
in case of strange attractors which exhibit multifractal structures, the convergence is 
slower. To perform the numerical simulations it has been used the L-moments proce- 
dure in order to overcome the difficulties of dealing with a singular continuous invariant 
measure. 



1 Introduction 

1.1 Classical Extreme Value Theory 

Extreme Value Theory (EVT), developed for the study of stochastical series of independent and 
identical distributed variables by Fisher and Tippett 19281 and formalized by Gnedenko 1943 , has 



been successfully applied to different scientific fields to understand and possibly forecast events that 
occur with very small probability but that can be extremely relevant from an economic or social 



point of view: extreme floods Gumbel 


1941 , Sveinsson and Boes, 2002 


, P. and Hense, 2007 , 


amounts of large insurance losses Brodin and Kluppelberg 2006|, |Cruz 


2002 ; extreme earth- 


quakes |Sornette et al. 1996|, |Cornell 


1968 , Burton 


1979|; meteorological and climate events 


Felici et al. 2007 , |Vitolo et al. 2009b , |Altmann et al. 


2006 , |Nicholis 1997 , Smith 1989 . An 



extensive review of the techniques and applications related to the EVT is presented in Ghil et al. 



2011 



Gnedenko [1943 studied the convergence of maxima of i.i.d. variables 

with cumulative distribution (cdf) F{x) of the form: 

F(x) = P{a rn {M m - b m ) < x} 



= max{Xo, X±, X m -i}. It may be rewritten 
Under general hypotesis on the nature of the 
parent distribution of data, Gnedenko 1943 show that the distribution of maxima, up to an affine 



Where a m and b m are normalizing sequences and M m 
as F(u m ) = P{M m < u m } where u m = x/a m + b m . 



change of variable, obeys to one of the following three laws: 
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Type 1 {Gumbel). 
Type 2 (Frechet). 

Type 3 (Weibull). 



E(x) = exp(— e x ), — oo < x < oo 



E(x) 



0, x < 



£(x) 



exp(— x £), for some £ > 0, x > 



exp(— (— x) for some £ > 0, a; < 
1, x > 



(2) 



(3) 



Let us define the right endpoint xf of a distribution function F(x) as: 

x,p = sup{x : F(x) < 1} 



(4) 



then, it is possible to compute normalizing sequences a m and b m using the following corollary 
of Gnedenko's theorem : 

Corollary (Gnedenko): The normalizing sequences a m and b m in the convergence of normalized 
maxima P{a m (M m — b m ) < %} — > F(x) may be taken (in order of increasing complexity) as: 

• Type 1: a m = [G(7 m )]"\ b m = -y m ; 

• Type 2: a m = 7m> b m = 0; 

• Type 3: a m = (x F - 7 m )~\ b m = x F ; 



where 



7 m = F^(l - l/m) = inf{x;F(x) > 1 - 1/m} 



G(t) 



XF l-F(u) 



du, t < xp 



(5) 



(6) 



In Faranda et al. 2011 we have shown that this approach is equivalent to fit unnormalized data 



directly to a single family of generalized distribution called GEV distribution with cdf: 



F G (x; a, £') = exp 



1 + ? 



x — fl 
a 



(7) 



which holds for l+£'(x— /J-)/cr > 0, using /* £ R (location parameter) and a > (scale parameter) 
as scaling constants in place of b m , and a m |Pickands III 1968] , in particular, in Faranda et al. |2011| 
we have shown that the following relations hold: 



/u = b % 



a 
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£' E R is the shape parameter also called the tail index: when £' — > 0, the distribution cor- 
responds to a Gumbel type ( Type 1 distribution). When the index is positive, it corresponds to 
a Frechet (Type 2 distribution); when the index is negative, it corresponds to a Weibull (Type 3 
distribution). 



To analyze the extreme value distribution in a series of data two main approaches can be applied: 
the Peak-over-threshold approach and the Block-Maxima approach. The former consists in look- 



ing at exceedance over high thresholds Todorovic and Zelenhasic, 1970] and a Generalized Pareto 



distribution is used for modeling data obtained as excesses over thresholds Smith 1984 , Davison 



19841, Davison and Smith 1990 



The so called block-maxima approach is widely used in climatological and financial applications 
since it represents a very natural way to look at extremes in fixed time intervals: it consists of 
dividing the data series of some observable into bins of equal length and selecting the maximum 



1992 , |Katz||1999|, |Katz et al.[ [2005 . 



(or the minimum) value in each of them Coles et al. 1999 , Felici et al. 2007|, Katz and Brown 



1.2 Extreme Value Theory for dynamical systems 



As far as the classical EVT is concerned, we should restrict our domain of investigation to the 
output of stochastic processes. Obviously, it is of crucial relevance for both mathematical reason 
and for devising a framework to be used in applications, to understand under which circumstances 
the time series of observables of deterministic dynamical system can be treated using EVT. 
Empirical studies show that in some cases a dynamical observable obeys to the extreme value 



statistics even if the convergence is highly dependent on the kind of observable we choose Vannitsem 


2007 1, Vitolo et al. 2009b , Vitolo et al. 2009a . For example, 


Balakrishnan et al. 


1995 and more 



recently Nicolis et al. 2006 and Haiman 2003 have shown that for regular orbits of dynamical 



systems we don't expect to find convergence to EV distribution. 
The first rigorous mathematical approach to extreme value theory in dynamical systems goes back 
to the pioneer paper by P. Collet in 20 01 |Collet 2001 . Important contributions have successively 



been given by Freitas and Freitas 2008 1 , Freitas et al. 2009].Freitas et al. 2010a and by Gupta et al. 



20091. The starting point of all these investigations was to associate to the stationary stochastic 
process given by the dynamical system, a new stationary independent sequence which enjoyed 
one of the classical three extreme value laws, and this laws could be pulled back to the original 
dynamical sequence. To be more precise we will consider a dynamical system (£l,B, v, f), where 
Q is the invariant set in some manifold, usually M d , B is the Borel cr-algebra, / : £1 — > Q, is a 
measurable map and v a probability /-invariant Borel measure. The stationary stochastic process 
given by the dynamical system will be of the form X m = g o f m 1 for any m £ N, where the 
observable g has values in MU±oo and achieves a global maximum at the point £ 6 $7. We therefore 
study the partial maximum M m = max{Xo, . . . ,X m _i}, in particular we look for normalising real 
sequences {a m },{b m },m G IR + for which i/{x;a m (M m — b m ) < t} = u{x\M m < u m } converge 
to a non-degenerate distribution function; here u m = + b m is such that mviX^ > u m ) — > t, 



for some positive r depending eventually on t: we defer to the book |Leadbetter et a~ 1983 



for a clear and complete picture of this approach. We will associate to our process a new i.i.d. 
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sequence Xq, ■ ■ ■ ,X m —i whose distribution is the same as that of Xq and with partial maximum: 
M m = max{Xo, • • • , X m -i}. Properly normalized the distribution of such a maximum converges 
to one of the three laws in equations T][3 and this is the interesting content of the Extreme Value 
Theory. Equations [T]{3] will be satisfied by our original process too, whenever we would be able to 
prove that 

lim u(M m < u m ) = lim v{M m < u m ) 

m— >oo m— >oo 

This can be achieved if one can prove two sufficient conditions called D2 and D' and which we briefly 
quote and explain in the footnote: these conditions basically require a sort of independence of the 
stochastic dynamical sequence in terms of uniform mixing condition on the distribution functions. 



In particular condition D2, introduced in its actual form by Freitas-Freitas Freitas and Freitas 2008 
could be checked directly by estimating the rate of decay of correlations for Holder observables [jj 
Another interesting issue of of the previous works was the choice of the observables g's: it is chosen 
as a function g(dist(x, £)) of the distance with respect to a given point £, with the aim that g 
achieves a global maximum at almost all points £ £ $7; for example g(x) = — logx. In particular the 
observable g was taken in one of three different classes gi,g2,g3, see Sect. 2 below, each one being 
again a function of the distance with respect to a given point £. The choice of these particular forms 
for the g's is just to fit with the necessary and sufficient condition on the tail of the distribution 
function F(u) = v{x;Xq < u}, in order to exist a non-degenerate limit distribution for the partial 
maxima I Freitas et al. 20091, iHolland et al. 2008 . We use here the fact that, thanks to conditions 



D2 and D' , the distributions of Xq rules out the distribution of our non-independent process X m 
as well. It is important to remind that the previous conditions will determine the exponent £ in the 
types 2 and 3 for E(x). 

Another major step in this field was achieved by establishing a connection between the extreme 



value laws and the statistics of first return and hitting times, see the papers by Freitas et al. 2009 



and Freitas et al. 2010b . They showed in particular that for dynamical systems preserving an 
absolutely continuous invariant measure or a singular continuous invariant measure v, the existence 
of an exponential hitting time statistics on balls around v almost any point £ implies the existence 
of extreme value laws for one of the observables of type gi, i = 1, 2, 3 described above. The converse 
is also true, namely if we have an extreme value law which applies to the observables of type 
gi,i = 1,2,3 achieving a maximum at £, then we have exponential hitting time statistics to balls 
with center £. Recently these results have been generalized to local returns around balls centered 



at periodic points Freitas et al. 2010a|. 



In the context of singular measures, the EVT has been developed in the recent paper by Freitas 



et al. |2010b . The main goal of their paper was to establish a connection with hitting and return 
time statistics; for that purpose they considered returns in balls and also into cylinders. We are 
particularly interested in their Theorem 1, about balls, since it covers the class of observables 
considered in this paper; in particular we use here one direction of the theorem which allows us 



1 We briefly state here the two conditions, we defer to the next section for more details about the quantities 
introduced. If X m , m > is our stochastic process, we can define Mjj = {Xj, Xj+i, ■ ■ ■ , Xj + {\ and we put Mo, m = 
Mm. The condition Z?2(Um) holds for the sequence X m if for any integer l,t,m we have |f(Xo > u m ,M t ,i < 
u m ) — v(Xq > Um)v(Mt,i < Um)\ < "/(m,t), where y(m, i) is non-increasing in t for each m and m-){m, t m ) — > as 
m — > 00 for some sequence t m — o(m), t m — s> 00. 

We say condition D'(u m ) holds for the sequence X m if lim;_ ,;r " *" vt m / ! ] 



limsup m mX)Ei v(X > u m ,Xj > u m ) = 0. 
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to get the extreme value distributions if the exponential return time statistics has been previously 
established for balls centered around almost all points and with respect to the probability invariant 
measure (in this manner we do not need to check conditions D2 and D'). 



1.3 This work 



In our previous work Faranda et al. [201 1 we have shown that there exists an algorithmic way to 
study EVT by using a block-maxima approach for dynamical systems which possess an absolutely 
continuous invariant measure and satisfy the mixing properties given by conditions D2 and D'. We 
have established the best conditions to observe convergence to the analytical results highlighting 
deviations from theoretical expected behavior depending on the number of maxima and number 
of block-observation. Furthermore, we have verified that the normalising process of variables can 
be applied a posteriori and a fit of unnormalised data produce a distribution that belongs to the 
Generalised Extreme Value (GEV) distributions family. 

In this work we test our algorithm for maps that do not have an absolutely continuous invariant 
measure. We remind that in the context of dynamical system, the invariant measure plays the role 
of the probability measure on the space of events; in this respect the general theory of extremes 
will continue to apply no matter such a probability is absolutely continuous or singular with respect 
to Lebesgue. The interesting point is that for the choice of observables we did (the functions c/j), 
the cumulative distribution function F will be related to the scaling of the measure v(B r (z)) of a 
ball B r (z) of radius r and centered at the point z, and such a scaling exponent turns out to be the 
Hausdorff dimension of the measure (also known as information dimension), when the point z is 
generically chosen. The experimental and accessible parameters of the GEV distributions will be 
explicitly expressed in terms of such a dimension. 

In order to get the values of £ and of a m and b m for finite m one should know how the measure of the 
ball B r (z) behaves as a function of r and of z and for measures which are not absolutely continuous. 
We notice that for absolutely continuous measure that approach works, at least in a few cases, and 
we quote our previous paper for that. Instead for singular continuous measures like those supported 
on Cantor sets, we are not aware of any analytic result allowing to get the few orders expansion 
of v(B T (z)). This will prevent us to compute rigorously the normalising constant a m for type 1 
observables g±; instead we will get the the limiting values of b m for type 1 and the limiting values 
of a m for type 2 and 3. Moreover we could not compute rigorously the exponent £. The values 
proposed for those non-rigorous constants are obtained by simply approximating v{B r (z)) with r . 
The agreement with the numerical simulations suggests that there were good choices and suggests 
also a direct proof of the EVT for our observables and with the normalising constants indicated by 
our heuristic analysis. 

As explained in the previous subsection we can either check the conditions D2 and D' or the 
existance of an exponential return time statistics. The latter is the case of iterated function systems 
considered in section 3.2: these are in fact given by expanding maps (since they verify the so-called 
open set condition) and the exponential return time statistics for balls could be proved, for instance, 



using the technique in Bessis et al. 1987 . It will be also the case for the Henon attractor with 
the parameters studied by Benedicks and Carlesson: for those parameters the attractor exists and 
carries an SRB measure; moreover very recently Chazottes and Collet established the Poissonian 
statistics for the number of visits in balls around generic points w.r.t. the SRB measure. Our 
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numerical computation will concern instead the usual Henon attractor. Finally, we will consider 



the Lozi attractor, and in this case we will quote the result by Gupta et al. 2009 , which proves 



the existence of the extreme value distributions for the observables constructed with the functions 
gi and for balls around almost any point w.r.t. the SRB measure. As a final remark, we stress that 
the results by Freitas-Freitas and Todd have been proved under the assumption that v{B r {C t )) is a 
continuous function of r: this is surely true for all the previous examples and such a condition will 
play a major role in our next considerations too. 

This work is organized as follows: in Section 2 we present the analytical results for the EVT in maps 
with singular measures deriving the asymptotic behavior of normalising sequences and parameters. 
In Section 3 we present the numerical procedure used for the statistical inference of the GEV 
distribution and the numerical experiments that we have carried out for both singular measures 
generated with Iterated Function Systems and maps with a less trivial measures such as the Baker 
transformation, Henon and Lozi maps. Eventually, in Section 4 we present our conclusion and 
proposal for future work. 



2 Extreme Value Theory for maps with singular measures 
2.1 Definitions and Remarks 

Let us consider a dynamical systems (f2, B, v, /), where Q is the invariant set in some manifold, 
usually M. d , B is the Borel cr-algebra, / : $7 — > Q, is a measurable map and v a probability /-invariant 
Borel measure. 

As we said in the Introduction and in order to adapt the extreme value theory to dynamical systems, 
we will consider the stationary stochastic process Ao,Xi, ... given by: 

X m (x) = <?(dist(/ m (x), 0) Vm g N (8) 

where 'dist' is a distance on the ambient space f2, £ is a given point and g is an observable 
function, and whose partial maximum is defined as: 

M m = max{X , X m _i} (9) 

The probability measure will be here the invariant measure v for the dynamical system. We will 
also suppose that our systems which verify the condition D2 and D' which will allow us to use the 
EVT for i.i.d. sequences. As we said above, we will use three types of observables gi,i = 1,2,3, 
suitable to obtain one of the three types of EV distribution for normalized maxima: 

51 (x) = -log(dist(x,C)) (10) 

52 (x) = dist(x, cr 1/o (11) 

g 3 (x) = C-dist(x,C) 1/a (12) 

where C is a constant and a > 6 M. 
These three type of functions are representative of broader classes which are defined, for instance, 



7 



throughout equations (1.11) to (1.13) in Freitas et al. [2009]; we now explain the reasons and the 
meaning of these choices. First of all these functions have in common the following properties: 
(i) they are defined on the positive semi-axis [0, oo] with values into RU {+oo}; (ii) is a global 
maximum, possibly equal to +oo; (iii) they are a strictly decreasing bijection in a neighborhood V 
of with image W. Then we consider three types of behavior which generalize the previous specific 
choices: 

Type 1: there is a strictly positive function p : W — > R such that Vy £ R we have 

gi 1 (s + yp(s)) _ 



lim 

s->Sl(0) 



Type 2: 52(0) = +00 and there exists f3 > such that Vy > we have 

lim ajHsy) = p 

Type 3: 53(0) = D < +00 and there exists 7 > such that Vy > we have 

limg3 -(D- 8y) = 

The Gnedenko corollary says that the different kinds of extreme value laws are determined by 
the distribution of 



F(u) = v(X Q < u) 



(13) 



and by the right endpoint of F, xp. 



We need to compute and to control the measure v{B r {C)) of a ball of radius r around the point 
£. At this regard we will invoke, and assume, the existence of the following limit 

logi/(B r (C)) 



lim 

r-s-0 



logr 



for £ chosen u — a.e. 



(14) 



Moreover we will assume that u(B r (Q) is a continuous function of r (see Freitas et al. |2009| for a 



discussion of this condition which shows that all the examples considered in our paper will fit it). 



When the limit (14) exists on a metric space equipped with the Borel c-algebra and a probability 
measure v, it gives the Hausdorff dimension of the measure or information dimension, defined as 
the infimum of the Hausdorff dimension taken over all the set of v measure 1 lYoung 1982 . This 



limit could be proved to exist for a large class of dynamical systems and whenever v is an invariant 
measure: let us indicate it with A without written explicitly its dependence on v. For example, 
for a very general class of one-dimensional maps with positive metric entropy, A is equal to the 
ratio between the metric entropy and the (positive) Lyapunov exponent of v |Ledrappier 1981 . For 



two dimensional hyperbolic diffeomorphisms, A is equal to the product of the metric entropy times 



the difference of the reciprocal of the positive and of the negative Lyapunov exponents [Young 



1982 1 . The information dimension is a lower bound of the Hausdorff dimension of the support of 



the measure v and it is an upper bound of the correlation dimension |Yakov 


1998|, |Hentschel and 


Procaccia 


1983 , |Grassberger 


1983|, |Bessis et al. 


1988|, |Bessis et al. 


1987 


, |Cutler and Dawson 



1989 
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2.2 Limiting behavior of the Extreme Value Theory parameters 

We summarize the three basic assumptions for the next considerations: 

• Assumption 1: our dynamical system verifies conditions (-D2) an d D' . 

• Assumption 2: the measure of a ball is a continuous function of the radius for almost all the 
center points; moreover such a measure has no atoms. 



Assumption 3: the limit (14) exists (and its value is called A) at almost all points £. 



Equipped with these conditions it is now possible to compute rigorously a few of the expected 
parameters for the three types of observables. 



Case 1: gi(x)= -log(dist(x,£)). Substituting equation M into equation 13 we obtain that: 



1 - F(u) = 1 - u(g(dist(x, ()) < u) 

= l-v(- log(dist(x, C)) < u) (15) 
= i/(dist(z,C) < e~ u ) = v{B e -u(Q) 

xf = sup{it; F(u) < 1} 

To use Gnedenko corollary it is necessary to calculate xp', in this case xf = +00 as we will explain 
in the proof below. 



According to Corollary 1.6.3 in Leadbetter et al. 1983 for type 1 a m = [G("f m )] 1 and b 



7 m = F 1 (1 — — ). We now show how to get the limiting value of 7 m ; a similar proof will hold for 
type II and III. 

Proposition 1. Let us suppose that our system verifies Assumptions 1,2,3 above and let us consider 
the observable g\; then: 

hm^ = A 

Proof 

By our choice of the observable we have: 1 — F(^ m ) = v{B e - lm (£)) = ^; since the measure is not 
atomic and it varies continuously with the radius, we have necessarily that 7 m —> 00 when m — > 00. 
Now we fix 5 > and small enough; there will be depending on 5 and on £, such that for any 
m > m S,C we have 

- hm < log v(B e - lm (C)) + A7 m < 5y m (16) 



Since logm — A7 m = — [log v{B e - lm (()) + A7 m ] and by using the bounds (16) we immediately have 

~hm < logm - A7 m < 5jm 

which proves the Proposition. 
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It should be clear that the previous proposition will not give us the value of j m and of b m , which 
is equal to j m for type I observables. We have instead a rigourous limiting behavior and we will 
pose in the following: 

7m = b m ~ — log m 

The values for finite m could be obtained if one would dispose of the functional dependence of 
v(B r (Q) on the radius r and the center £: this has been achieved for non-trivial absolutely continuous 



invariant measure in our previous paper |Faranda et al. 2011] . The same reason prevent us to get 
a rigorous limiting behavior for a m = [G(7 m )] _1 . The only rigorous statement we can do is that 
G{^ m ) = o(7m); this follows by adapting the previous proof of the proposition to another result 
(see |Leadbetter et al. 1983] ) which says that for type I observables one has lim m _ 5 . 00 n(l — F{j m + 



x G(jm)}) = e x , for all real x: choosing x = 1 gives us the previous domination result. In the 
following and again for numerical purposes we will take 



This follows easily by replacing in formula ^ u(B r (()) ~ r A for r small. We finish this part by 
stressing that for our observable we expect £ = £' = 0. 

Case 2: g2(x)=dist(x,(")~ 1 / a . In this case we have 

1 - F(u) = 1 - i/(dist(s, Cr 1/Q < u) 

= 1 - i/(dist(x, C) > u~ a ) (17) 

and xf = +oo. Since b m = we have only to compute a m which is the reciprocal of 7 m which is in 
turn defined by 7 m = i ?_1 (l — 1/m). By adapting Proposition [l] we immediately get that 

logm 
lim = aA 

m^co log 7 m 

which we allow us to use the approximation a m ~ — t— The exponent ^ for Type II observables is 
given by the following limit (see |Leadbetter et al. 1983 , Th. 1.6.2) 



lim (1 - F(tx))/(l - F(t)) = x^, £ > 0, x > 



t— >co 



The crude approximation u(B r (()) ~ r A for r small, will give immediately that £ ~ aA and this 
value will appear in the exponent £' = 1/^ in the distribution function given by the GEV. 

Case 3: g 3 (x)=C-dist(x,C) 1/a . We have first of all: 

1 _ F{u) = i- u ( C - dist(x, C) 1/Q < u) 
= u(B {c _ u)a (0) 
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In this case xp = C < oo and a m = (C — 7 m ) ; b m = C. The previous proposition immediately 



log m 



A which gives the asymptotic scaling 7 m ~ C- 



shows that liim^ _ alog(c _ 7m) 

b m = C. Finally the exponent £ is given again by Th. 1.6.2 in Leadbetter et al. [1983 by the formula 

lim(l - F(C - hx))/(l - F(C -h)) = x i , £ > 0, x > 
h-+0 



which with our usual approximation furnishes £ ~ a A. 



3 Numerical Experiments 



3.1 Procedure for statistical inference 

For the numerical experiments we have used a wide class of maps that have singular measure, also 
considering the case of strange attractors such as the ones observed by iterating Lozi Map or Henon 

|201l] 



map. The algorithm used is the same described in Faranda et al 



for each map we run 

a long simulation up to k iterations starting from a given initial condition. From the trajectory 
we compute the sequence of observables gi,g2, 93 dividing it into n bins each containing m = k/n 
observations and eventually obtaining the empirical cdf of maxima. 



In Faranda et al. |201l| we have used a Maximum Likelihood Estimation (MLE) procedure working 
both on pdf and cdf (cumulative distribution function), since our distributions were absolutely 
continuous and the minimization procedure was well defined. In this case, we don't have anymore 
the pdf and consequently the fitting procedure via MLE could give us wrong results. To avoid these 
problems we have used an L-moments estimation as detailed in Hosking 1990 . This procedure 



is completely discrete and can be used both for absolutely continuous or singular continuous cdf. 
The L-moments are summary statistics for probability distributions and data samples. They are 
analogous to ordinary moments which meant that they provide measures of location, dispersion, 
skewness, kurtosis, but are computed from linear combinations of the data values, arranged in 
increasing order (hence the prefix L). Asymptotic approximations to sampling distributions are 
better for L-moments than for ordinary moments |Hosking |1990|, Figure 4]. The relationship 



between the moments and the parameters of the GEV distribution are described in Hosking (1990 



while the 95% confidence intervals has been derived using a bootstrap procedure. As comparison 
we have checked that the results presented in Faranda et al. [201 1 are comparable with L-moments 
methods. We have found that both the methods give similar results even if L-moments has an 
uncertainty on the estimation of parameters generally slightly bigger. 

The empirical cdf contains plateaux which correspond to non accessible distances in correspondence 
of the holes of the Cantor set. The discrete nature of L-moments allows to overcome difficulties 
that may arise in singular continuous cdf: the normalization procedure carried out with this method 
consists in dividing each quantity computed via L-moments by a function of the total number of 
data will prevent us from obtain a unnormalised distribution. The last issue we want to address is 
the choice of a suitable model for our data: in principle, using L-moments procedure we can fit the 
data to any kind of known cdf. To validate the use of the GEV model we proceeded in the following 
way: 

-A priori the choice of a GEV model arises naturally if the assumptions presented in the section 
3 are satisfied. In this set up we can directly compute the parameters of GEV distribution using 
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L-moments as described in Hosking 1990 . 



-A posteriori we can verify the goodness of fit to GEV family if we apply some parametric or non- 
parametric tests commonly used in statistical inference procedures. For this purpose we have fitted 
our experimental data to a wide class of well known continuous distributions. Using Kolmogorov 



Smirnov test (see Lilliefors 1967 for a description of the test) we have measured the deviation 



between the empirical cdf and the fitted cdf, finding that using the GEV distribution we effectively 
achieve a minimization of the deviation parameter. 

We summarize below the results we expect from numerical experiments in respect to n according 
to the conjecture described in the previous section. Since we keep the length of the series k = n - m 
fixed, the following relationships can be obtained simply replacing m = k/n in the equations derived 
in the previous section. 
For gi type observable: 



a 



1 

A ^ A 

For 52 type observable, we can either choose b 



]n(k/n) £' = 



(19) 

where c G 1 is positive 

constant, as detailed in Beirlant |2004|. A priori, we do not know which asymptotic sequences will 



or b r , 



c ■ m 



correspond to the parameters \x in the experimental set up. The experimental procedure we use 
automatically select b m = c • , therefore the following results are presented taking into account 
this asymptotic sequence: 



a ~ n 



-l/(aA) 



\i ~ n 



-l/(oA) 



i 

aA 



For 53 type observable: 



a ~ n 



l/(aA) 



H = C £' 



1 
aA 



(20) 



(21) 



3.2 IFS for Cantor Sets 

A Cantor set can be obtained as an attractor of some Iterated Function Systems (IFS). An IFS is a 
finite family of contractive maps {/i, fa, f s } acting on a compact normed space f2 with norm | • | 
and possessing a unique compact limit set (the attractor) Kg!! which is non-empty and invariant 
by the IFS, namely: 



K 



U /<(*)• 



We will put a few restrictions on the IFS in order to see it as the inverse of a genuine dynamical 
system; we will explain in a moment why this change of perspective will help us to compute observ- 



ables on fractal sets. We defer to the fundamental paper by Barnsley and Demko 1985] , for all the 
material we are going to use. 

First of all we will consider the /j as strict contractions, namely there will be a number < A < 1 

such that for all i = 1, • • • , s we have \ fi(x) — fi(y)\ < \\x — y\, for all x, y G X. 

Then we will suppose that each /, is one-to-one on the attractor K and moreover Vz = 1, • • ■ , s 
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we have fi(K) n fj(K) = 0, i 7^ j (open set condition). This will allow us to define a measurable 
map T : K — > K by T(x) = ff 1 (x) for x £ fi(K): the attractor if will be the invariant set for 
the transformation T which will play therefore the role of a usual dynamical system. A complete 
statistical description of a dynamical system is given by endowing it with an invariant probability 
measure; in particular we ask that this measure by ergodic if we want to compute the maxima of 
the sequence of events constructed with the observables g±\ we remind that these events are nothing 
but the evaluations along the forward orbit of an initial point chosen according to the measure. If, 
as always happens, the attractor K has a fractal structure and zero Lebesgue measure, we could not 
get such an initial point for numerical purposes. We overcome this situation for attractors of global 
diffeomorphisms by taking the initial point in the basin of attraction and by iterating it: the orbit 
will be distributed according to the SRB measure (we will return later on this measure). For our 
actual attractors generated by non- invert ible maps, the iteration of any point in the complement of 
the attractor will push the point far from it: it would be better to call repellers our invariant sets 
instead of attractors. We have therefore to proceed in a different manner. The measures supported 
on the attractor K will give the solution. First of all let us associated to each map fi a positive 
weight pi in such a way that J2i=iPi = 1- Then it is possible to prove the existence of a unique 
measure v (called balanced) which enjoy the following properties: 

• The measure v is supported on the attractor K and it will be invariant for the map T associated 
to our IFS (see above). 

• For any measurable set B in X we have 



• Let us put (Sg)(x) = Yli=iPi9(fi( x ))> f° r a continuous functions g on £1 and for any point 
xell, then we have 



This last item is very important for us; first of all it holds also for the characteristic function of a 
set provided the boundary of this set has v measure zero; therefore it is a sort of ergodic theorem 
because it states that the backward orbit constructed by applying to any point in X the maps fi 
with weights pi will distribute on the attractor K as the forward orbit (namely the orbit generated 
by the transformation T associated to the IFS) of a point y G K and chosen almost everywhere 
according to v. 
Let us give our first example. 

3.2.1 Uniform weights and the Sierpisnkij triangle 

We consider the middle one third Cantor set that is the attractor of the IFS {/1, j{\ defined as: 



s 



(B) = ^v{f-\B)) 



i=i 




fi(x) = x/3 with weight p\ 
H{ x ) = (x + 2)/3 with weight pi 



(22) 
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where x E [0, 1] and we set p\ = p 2 = 1/2 so that, at each time step, we have the same probability 
to iterate fi(x) or ^(x): 

Equivalently the previous IFS can be written as: 



x t+ i = (x t + b)/3 (23) 
where, at each time step, we extract randomly with equal probability b to be or 2. 

We will consider also the so called Sierpisnkij triangle, defined by 



x t +i = (x t + u p ,i)/2 

yt+i = {yt + v P)2 )/2 



(24) 



We extract randomly at each time step and with equal probability the number p to be the integer 



1,2 or 3. Then we iterate the map 24 substituting the elements v p> i and v Pt 2 of the following matrix: 



1 
-1 
1 

For these attracting sets the information dimensions are well known, they are A 
for the Cantor set and A = log(3)/log(2) for the Sierpinskij triangle |Sprott 2003 



log(2)/log(3) 



In the following experiments In order to choose the centers £ of our balls we proceed by using again 
the backward preimages of any point in f2, namely we take a point x 6 f2 and we consider £ as one 
of the preimages f~ t (x) with t much larger than the sequence of observed events; by what we said 
above, that preimage will be closer and closer to the invariant Cantor set and also it will approaches 
a generic point with respect to the balanced measure v. 

First of all we have analised the empirical cdf F(u) of the extrema for g\ observable. An example 
is shown for the IFS in eq. 22 in Figure [T] The histogram is obtained iterating the map in equation 



22 for 5-10 iterations, Q ~ 0.775, a = 4, C = 10. Once computed the functions g\ the series of 



maxima for each observable is computed taking each of them in bins containing 5000 values of g\ 
for a total of 1000 maxima. As claimed in the previous section, the cdf is a singular continuous 
function and this is due to the structure of the Cantor set. The results are similar for the other 
observables and other initial conditions. 
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Figure 1: Empirical (blue) and fitted (red) cdf for IFS that generates a Cantor Set, £ = 0.775, g\ 
observable function. 



To check that effectively the parameters of GEV distribution obtained by L-moments estimation 
are related to the fractal dimension of the attracting Cantor Set and the Sierpinskij triangle, we 
have considered an ensemble of 10 4 different realizations of the eq. 22 and eq. 24, starting from the 



same initial conditions. To check the behavior we have varied n and m keeping fixed the length of 
the series k = 10 7 . 



In 



Faranda et al. 2011 we have shown that a good convergence is observed when 



n, m > 1000, therefore we will make all the considerations for (n, m) pairs that satisfy this condition. 



In figures |2j|4| the results of the computation for the IFS that generates the Cantor Set ( plots on 
the left) and the Sierpisnkij triangle (plots on the right) are presented. In all the cases considered 
the behavior well reproduce the theoretical expected trend described in equations 1~9]|2~T The initial 
condition here shown is £ = 1/3 for Cantor Set and Q — (0.02,0.40) for Sierpinskij triangle, but 
similar results hold for different initial conditions if chosen on the attractive sets. The black line 
is the mean value over different realizations of the map, while the black dotted lines represent one 
standard deviation. 



For g\ observable function, according to equation 19 



we expect to find £ = in both cases 
and this is verified by experimental data shown in Figure]^,). For the scale parameter a similar 
agreement is achieved in respect to the theoretical parameters A = a ^ gi ^ = l°g(2)/ log(3) — 0.6309 
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A = l/(a|«|) 


Cantor 


Sierpinskij 


Theoretical 


log(2)/log(3) ~ 0.6309 


log(3)/log(2) ~ 1.5850 




0.636 ±0.006 


\k\ = 1.592 ±0.007 


a {92) 


0.634 ±0.007 


|k| = 1.56 ±0.02 


c(93) 


0.64 ±0.01 


\k\ = 1.62 ±0.01 



Table 1 : Information dimension A computed taking the logarithm of equations 20]j2T and computing 
the angular coefficient k of a linear fit of data; IFS with uniform weights and Sierpinskij triangle, 
a = 4, C = 10. 



for the Cantor Set and A 



a ( gi ) = log(3)/ log(2) ~ 1.5850 for Sierpinskij shown in figure 2p) with 
a green line. Eventually, the location parameter \i shows a logarithm decay with n as expected from 

The linear 
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equation 

fit angular coefficient k of equation 19 satisfies A 
using this relation is A 



A linear fit of fi in respect to log(n) is shown with a red line in figure [2j: 

1/k and the dimension computed from data 
1.59 ± 0.01 in the Sierpinskij triangle. 



0.64 ± 0.01 for Cantor and A 
The agreement between the conjecture and the results are confirmed even for 52 type and 53 type 
observable functions shown in figures [3] and [4] respectively. In this case we have experienced some 
problems in the convergence for the Cantor map has using a = 2 and a = 3. The problem is 
possibly due to the fact that L-moments method works better if £' £ [—0.5, 0.5] while for a < 3 the 
shape parameter > 0.5 for the Cantor map. For this reason results are shown using a = 4 both 
for Sierpinskij triangle and Cantor IFS, for all the experiments presented the constant value in 53 
will be C = 10. For both observables 52)53 there is strong agreement between the experimental 
and theoretical £' values. In figures]^), [3^), |4j)) a log-log scale is used to highlight the behavior 
described by eq. 20 and 21 We can check again the value of the dimension using the angular 
coefficient k which satisfies A = l/(a|/-t|). In Table 1 we compare these results with the theoretical 
values. The error is here represented as one standard deviation of the ensemble of realizations and 
we find a good agreement between theoretical and experimental parameters within two standard 
deviations. Eventually, computing 53 observable function we expect to find a constant value for 
fj,{g$) = C = 10 while a{g^) has to grow with a power law in respect to n as expected comparing 



with equation 21 



Other tests have been done computing the statistics using parameter a 
observables. Also in this cases, no deviation from the behavior described in eq. 
found. 



5,6,7,8 for g 2 and g 3 
19]|2T] has been 
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Figure 2: g\ observable, a) £' VS log 10 (n); b) a VS log 10 (n); c) \x VS log(n). Cantor set, Right: 
Sierpinskij triangle. Dotted lines represent one standard deviation, red lines represent a linear fit, 
green lines are theoretical values. 17 



(a) 




Iog10(n) Iog10(n) 



Figure 3: g 2 observable a) £' VS log 10 (n); b) log 10 (<r) VS log 10 (n); c) log 10 (/x) VS log 10 (n). Left: 
Cantor set, Right: Sierpinskij triangle. Dotted lines represent one standard deviation, red lines 
represent a linear fit, green lines are theoretical iSlues. 




Figure 4: g 3 observable, a) £' VS log 10 (n); b) log 10 (<r) VS log 10 (n); c) log 10 (» VS log 10 (n). 
Cantor set, Right: Sierpinskij triangle. Dotted lines represent one standard deviation, red lines 
represent a linear fit, green lines are theoretical ifilues. 



3.2.2 IFS with non-uniform weights 

Let us now consider the case of an IFS with different weights: 



f k (x) = a k + X k x £€[0,1] fe = l,2, ...,s 



(25) 



and each fa is iterated with (different) probability lOj. 

In this case it is possible to compute the information dimension as the ratio between the metric 
entropy and the Lyapunov exponent of the associated balanced measure Barnsley 2000 : we get 
the following expression: 



A 



w\ \ogwl + ... + w s log w s 



wi logAi + ... + w s logA s 
In this analysis we have considered the following IFS: 




x/3 

(x + 2)/3 



with weight w 
with weight 1 - 



w 



(26) 



(27) 



and we have changed the weight w between 0.35 and 0.65 with 0.01 step. For w = 0.5 we 
obtain the same results shown in the previous section, while for different weights we can check the 



expression 26 



The presence of different weights makes the convergence process sensible to the choice of the 
sample point £ where our observable reaches its maximum. For that reason we took several different 
values of £ in order to obtain a reliable estimations of the information that should be obtained, in 
this case, as an average property. 

In Figure [5] we present the dimension A computed using relationship 19 21 In particular we 
can compute the dimension from eq. 19 as: 



AH 51 )) 



We can infer dimension also from eq. 20 eq. [21] as: 



A(e'(^)) 



i 



> 



2,3 



(28) 



(29) 



< eta 

and in all expressions above the brackets < . > indicate an average on different sample points 
(. For the rest of the numerical computations we set a = 5. The parameters have been computed 
using 1000 different initial conditions on the support of the attractor, and for 30 realizations of each 
sample point £,the block-maxima approach is here used with n = m = 1000. The error bar are 
computed using the standard error propagation rules. 

The agreement between the theoretical dimension and the experimental data is evident for 
all the weights and for all the observable considered. The uncertainty increases when w is much 
different from 0.5. This is due to the fact that as soon as we change the weight to be different from 
0.5 the parameters spread increase to take in account the local properties of the attractor. The 
best agreement and less uncertainty is achieved considering the dimension as computed from cr(gi) 
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observable. This is possibly due to the slower convergence for gi or g% observables to the respective 
theoretical distributions: in g\ we modulate the distances with a logarithm function while in c/2 and 
in 53 power laws are used. Nevertheless, all the data show the right trend. 



0.68 - 



0.66 



0.64 - 



0.62 




0.58 



0.56 - 



0.54 



0.52 



Figure 5: Different estimation of A dimension obtained using Extreme Value distribution for the 



IFS in equation 27 



3.3 Non-trivial singular measures 

In the previous subsection we have analysed the relatively simple cases of Cantor sets generated 
with IFS. In order to provide further support to our conjectures, we now present some application 
of our theory to the output of dynamical systems possessing a less trivial singular measures. We 
consider three relevant examples of two dimensional maps. 
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The Baker map 

The Baker map is denned as follows: 



x t +i = 



\ j a x t mod 1 if y t < a 

\l/2 + -f b x t modi if y t > a 



mod 1 if y t < a 



I mod 1 if y t > a 



we consider the classical value for the parameter: a = 1/3, 7 = 1/5 and % = 1/4. 

Rigorous analytical results are available for the estimation of the information dimension |Kaplan 



and Yorke 1979 . For our parameter values, the analytical expected value is D ~ 1.4357. 

We have performed the same analysis detailed in section 4.1, but with a difference. This map 
is invertible and its invariant set is an attractor given by the cartesian product of a segment, along 
the y-axis, and a one dimensional Cantor set along the x-axis. The system possess an invariant 
SRB measure, which can practically be constructed by taking ergodic sums for any point sitting on 
the basin of attraction. In order to compute the center of the balls on the attractor, we proceed in 
a similar manner as for repellers (see above), namely we take any point x in the basin of attraction 
and we iterate it t times with t much bigger than the sequence of observed events. Then we take 
£ as the point /*(#): it will be closer and closer to the attractor and distributed according to the 
SRB measure. In our set up: a = 4 for gi and g%, C = 10. The results are shown in figures [6]|8] the 
black continuous lines will represent the parameter average over different initial conditions and the 
black dotted lines the standard deviation of the distribution of the estimated parameters. 
The expected theoretical £' values are within one standard deviation of the results of the fit for all 
three observable. The agreement seems to be better when we increase n even if this correspond to a 
decrease of m in our set up. This behavior is quite interesting since it seems that we obtain a much 
better convergence to theoretical values if n ~ 10 4 , while in all other example there is no such a dif- 
ference between n = 1000 and m = 10000. A similar consideration can be made for cr(gi) shown in 
Figure]^)) that approaches the theoretical values value for bigger n values. The angular coefficient 
of the linear fit for fi(gi) shown in the semilog plot in Figure ^p) allow us to estimate the dimension 



A = 1/\k\ = 1.48±0.03 that is consistent with the theoretical values within two standard deviations . 

Log-log plots of the parameters against n are shown in figures [Tja) , |7p) and ^p) , and the value of 
the dimension A computed using the angular coefficients are reported in Table 2. The agreement 
with expected value is good enough for all the parameters and better for ^{g^j- Eventually, Figure 
[8^) shows that ix{gz) approaches C = 10. 
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Figure 6: g\ observable, a) £' VS log 10 (n); b) a VS log 10 (n); c) fi VS log(n). Baker map. Dotted 
lines represent one standard deviation, red lines represent a linear fit, green lines are theoretical 
values. 23 




Figure 7: <?2 observable a) £' VS log 10 (n); b) log 10 (cr) VS log 10 (n); c) log 10 (/x) VS log 10 (n). Baker 
map. Dotted lines represent one standard deviation,, red lines represent a linear fit, green lines are 
theoretical values. 24 
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Figure 8: # 3 observable, a) £' VS log 10 (ra); b) log 10 (a) VS log 10 (n); c) log 10 (/x) VS log 10 (n). Baker 
map. Dotted lines represent one standard deviation, red lines represent a linear fit, green lines are 
theoretical values. 25 



The Henon and Lozi maps 

The Henon map is denned as: 



x t +i 
Vt+i 



bxt 



while in the Lozi map ax\ is substituted with a\xt\ 



xt+l = Vt + 1 
2/t+i = ^ 

We consider the classical set of parameter a 
- 0.5 for the Lozi map. 



a\xt\ 



(31) 



(32) 



1.4, 6 = 0.3 for the Henon map and a = 1.7 and 



Young 1 1985 proved the existence of the SRB measure for the Lozi map, whereas for the Henon 



map no such rigorous proof exists, even if convincing numerical results suggest its existence |Badii 



and Politi, 1987]. Note that Benedicks and Carleson 1991 proved the existence of an SRB measure 



for the Henon map with a different set of parameters. Using the classical Young results which makes 
use of the Lyapunov exponents, we obtain an exact result for A for the Lozi attractor: 

A ~ 1.40419 

Instead, in the case of the Henon attractor, we consider the numerical estimate provided by 
Grassberger |1983 : 



A = 1.25826 ±0.00006 

As in the previous cases the GEV distribution is computed with L-moments methods varying 
n and m and averaging the distribution parameters over 1000 different sample points chosen as 
described before for the Baker map. Results are presented in figures 9]|TT the plots on the left-hand 
side refer to the Henon map, while on the right-hand side the results refer to the Lozi map. 
When considering the numerical results are in agreement with the theoretical estimates. Never- 
theless, the parameters distribution have a rather range spread which indicates a slower convergence 
towards the expected values in respect to what is observed for the IFS case. The experimental val- 
ues of cr(gi) , shown in Figure [9]} , approach the theoretical values shown by a green line. The 
angular coefficient computed from the semilog plot of n{g\) represented in Figure [9^) gives us an 
estimate of the dimension A = 1/|k| = 1.234±0.015 for Henon and A = l/lfd = 1.40±0.01 for Lozi. 



The other angular coefficients related to g2 and g^ observables for the plots shown in figures |l Op), 
10;), and |llb) are presented in Table 2. Within 2 standard deviation they are comparable with 
the theoretical ones and the best agreement is achieved considering /i(<?2)- The constant C = 10 is 



approached quite well (Figure 11 z)). 
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A = l/(a|«|) 


Baker 


Henon 


Lozi 


Theor. 


1.4357 


1.2582 


1.4042 




1.47 ±0.02 


1.238 ±0.009 


1.396 ± 0.008 




1.39 ±0.04 


1.35 ±0.07 


1.38 ±0.02 


0(93) 


1.56 ±0.08 


1.15 ±0.07 


1.42 ±0.01 



Table 2: Information dimension A computed taking the logarithm of equations 20 21 and computing 
the angular coefficient k of a linear fit of data; for Baker, Henon and Lozi maps. 



The slower convergence for these maps may be related to the difficulties experienced computing 



the dimension with all box-counting methods, as shown in Grassberger 1983 .Badii and Politi 1987 



In that case it has been proven that the number of points that are required to cover a fixed fraction 
of the attractor support diverges faster than the number of boxes itself for this kind of non uniform 
attractor. In our case the situation is similar since we consider balls around the initial condition £. 
As pointed out, the best result for the dimension is achieved using the parameters provided by g\ 
observable since the logarithm modulation of the distance exalts proper extrema while weights less 
possible outliers. 



27 




Figure 9: g\ observable, a) £' VS log 10 (n); b) a VS log 10 (n); c) \i VS log(n). Left: Henon map, 
Right: Lozi map. Dotted lines represent one standard deviation, red lines represent a linear fit, 
green lines are theoretical values. 28 




Figure 10: 52 observable a) £' VS log 10 (n); b) log 10 (<r) VS log 10 (n); fc) log 10 (//) VS log 10 (n). Left: 
Henon map, Right: Lozi map. Dotted lines represent one standard deviation,, red lines represent a 
linear fit, green lines are theoretical values. 29 




Figure 11: g 3 observable, a) £' VS log 10 (n); b) log 10 (cr) VS log 10 (n); c) log 10 (/i) VS log 10 (n). Left: 
Henon map, Right: Lozi map. Dotted lines represent one standard deviation, red lines represent a 
linear fit, green lines are theoretical values. 30 



4 Final Remarks 



Extreme Value Theory is attracting a lot of interest both in terms of extending pure mathematical 
results and in terms of applications to many fields of social and natural science. As an example, in 
geophysical applications is crucial to have a tool to understand and forecast climatic extrema and 
events such as strong earthquakes and floods. 

Whereas the classical Extreme Value Theory deals with stochastic processes, many applications 
demanded to understand whereas it could rigorously be used to study the outputs of deterministic 
dynamical systems. The mathematical models used to study them present a rich structure and their 
attracting sets are very often strange attractors. In such sense is extremely important to develop 
an extreme value theory for dynamical system with singular measures. Recently, The existence of 
extreme value laws for dynamical systems preserving an absolutely continuous invariant measure or 
a singular continuous invariant measure has been proven if strong mixing properties or exponential 
hitting time statistics on balls are satisfied. 

In this work we have extended the results presented in Faranda et al. 2011 to the case of dynamical 
systems with singular measures. Our main results is that there exist an extreme value distribution 
for this kind of systems that is related to the GEV distribution when observable functions of the 
distance between the iterated orbit and the initial conditions are chosen. The three extreme value 
type for the limit distribution laws for maxima and the generalized distribution family (GEV) are 
absolutely continuous function. We will recover the GEV using histograms on the frequency of 
maxima; in this way the cumulative distribution function which we got from such an histogram 
will have plateaux just in correspondence of the holes of the Cantor set, whenever this one is the 
invariant set. This could be easily explained by the very nature of our observables which measure 
the distance with respect to a given point: there will be distances which are not allowed when 
such distances are computed from points in the holes. It should be stressed that such a cumulative 
distribution function, which is a sort of devil staircase and therefore is a singular continuous function, 
in any way could converge to a GEV distribution. The latter as the three type extreme values 
laws are normalized laws which must be adjusted in order to give a non-degenerate asymptotic 
distribution. The strength of our approach, as we said above, relies in the possibility to infer the 
nature and the value of such normalizing constants by a fitting procedure on the unnormalized data, 
the histograms. This worked very well for probability measure which were absolutely continuous. 
We will see that it works also for singular measures (and the normalizing constants will be related 
to the information dimension), provided we remind that this time the fitting procedure will contain 
a sort of extrapolation to smooth out the gaps of the Cantor sets. 

It is interesting to observe that on Cantor sets the notion of generic point is not so obvious as for 
smooth manifolds which support Lebesgue measure: in this case in fact one could suppose that 
each point accessible for numerical iterations is generic with respect to an invariant measure which 
is in turn absolutely continuous. This notion of genericity is restored on attr actor by considering 
the SRB measure. Instead for Iterated Function Systems we can dispose of uncountably many 
measures, but we have a precise manner to identify them and this will be reflected in the different 
dimensions produced by the numerical computation of the parameters of the GEV. The possibility 
to discriminate among different singular measures having all the same topological support is another 
indication of the validity and of the efficiency of our approach. 

We have also shown that the parameters of the distribution are intimately related to the information 
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dimension of the invariant set. We have tested our conjecture with numerical experiments on 
different low dimensional maps such as : the middle third Cantor set, the Sierpinskij triangle), 
Iterated Function System (IFS) with non-uniform weights, strange attractors such as Lozi and 
Henon. In all cases considered there is agreement between the theoretical parameters and the 
experimental ones. The extimates of A are in agreement with the theoretical values in all cases 
considered. It is interesting to observe that the algorithm described with the selection of maxima 
acts like a magnifying glass on the neighborhood of the initial condition. In this way we have both 
a powerful tool to study and highlight the fine structure of the attractor, but, on the other hand we 
can obtain global properties averaging on different initial conditions. Even if we are dealing with 
very simple maps for which many properties are known it is clear from numerical experiments that 
is not so obvious to observe a good convergence to the GEV distribution. Even if we are able to 
compute very large statistics and the results are consistent with theoretical values, the error range 
is wide if compared to the experiment for maps with a.c.i.m. measures that we have carried out in 



Faranda et al. 2011 . This should be taken in consideration each time this statistics is applied in a 
predictive way to more complicated systems. 

In the case of an experimental temporal series, for which the underlying dynamics is unknown, a 
classical problem is to obtain the dimensionality of the attractor of the dynamical systems which 
generated it. This can be achieved through the so called Ruelle-Takens delay embedding where, 
starting from the time series of an observable 0(n), we can construct the multivariate vectors in a 
A-dimensional space: 

4>{N) = [0(n),0(n + 1), 0(n + d- 1)] 



and study the geometrical properties using the Recurrence Qualification Analysis Marwan et al. 



20071. The minimum value of A needed to reconstruct the actual dimension d* is given by [2d] + 2. 



Using the procedure described in this paper we could find an estimate for d* thus determining the 
minimum value of A to be used in the Recurrence Quantification Analysis. We will test this strategy 
in a subsequent publication. 



The theory and the algorithm presented in this work and in Faranda et al. 2011 allow to study in 
detail the recurrence of an orbit around a point: this is due to the particular choice of the observables 
that require to compute distances between initial and future states of the system. Understanding 
the behavior of a dynamical system in a neighborhood of a particular initial condition is of great 
interest in many applications. As an example, in weather forecast and climate it is important to 
study the recurrence of patterns (the so called analogues). In principle, applying the extreme value 
statistics to the output of meteorological models, will make possible to infer dynamical properties 
related to the closest return towards a certain weather pattern. EVT will give information not only 
about the probability distribution of the extrema but also about the scaling of the measure of a 
ball centered on the chosen initial condition providing an insight to the dynamical structure of the 
attractor. 
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